/*
 * intepolate.cpp
 *
 *  Created on: May 30, 2010
 *      Author: snirgaz
 */

#include "Interpolate.h"
//#define FORCE

Interpolate::Interpolate() {
	String ds_name("Dataset1");
	dr_.resize(QmcParams::NumP, QmcParams::M, QmcParams::dim);
	for (int d = 0; d < QmcParams::dim; d++) {
		resFProposed_[d].resize((QmcParams::mbar + 3) * QmcParams::NumP);
		resFOriginal_[d].resize((QmcParams::mbar + 3) * QmcParams::NumP);
	}
	resV_.resize((QmcParams::mbar + 3) * QmcParams::NumP);
	resF_.resize((QmcParams::mbar + 3) * QmcParams::NumP);
#ifdef POT_LOG
	LogArg_ = ippsMalloc_64f((QmcParams::mbar + 3) * QmcParams::NumP);
#endif
#ifdef POT_AZIZ
	drSqr_.resize((QmcParams::mbar + 3) * QmcParams::NumP);
#endif
#ifdef EXT_POT
	QmcParams::gridResExtPot = int(double(QmcParams::NumP) * QmcParams::disXi);
	dxExtPot_ = 1 / (double(QmcParams::gridResExtPot - 1));
	dResExtPot_ = (double(QmcParams::gridResExtPot - 1));
	extPotFOriginal_.resize(QmcParams::NumP);
	extPotFProposed_.resize(QmcParams::NumP);
#ifdef STEPS
	vgridExtPot_ = new double[QmcParams::gridResExtPot
	* QmcParams::gridResExtPot];
	unsigned int seed = QmcParams::seed;
	ippsRandUniform_Direct_64f(vgridExtPot_,
			QmcParams::gridResExtPot * QmcParams::gridResExtPot,
			-QmcParams::disStrength, QmcParams::disStrength, &seed);
#endif
#if not(defined(LATTICE) or defined(STEPS))
	vgridExtPot_ = new double[QmcParams::gridResExtPot
	* QmcParams::gridResExtPot];for
	( int d = 0; d < SIM_DIM; d++)
	fgridExtPot_.push_back(new double[QmcParams::gridResExtPot * QmcParams::gridResExtPot]);
	hdf5_rw::read_double(QmcParams::GridFileVExtPot, ds_name, vgridExtPot_);
	hdf5_rw::read_double(QmcParams::GridFileFXExtPot, ds_name, fgridExtPot_[0]);
	hdf5_rw::read_double(QmcParams::GridFileFYExtPot, ds_name, fgridExtPot_[1]);
#endif
#endif
}
Interpolate::~Interpolate() {
#ifdef POT_LOG
	for (int d = 0; d < SIM_DIM; d++) {
		ippsFree(resFProposed_[d]);
		ippsFree(resFOriginal_[d]);
	}
	ippsFree(resV_);
	ippsFree(LogArg_);
#endif
#ifdef STEPS or EXT_POT and not defined (LATTICE)
	delete[] vgridExtPot_;
#endif
#ifdef not defined (STEPS)
	for (int d = 0; d < SIM_DIM; d++) {
		delete[] fgridExtPot_[d];
	}
#endif
#ifdef POT_AZIZ
	//ippsFree(resV_);
	//ippsFree(resF_);
	//ippsFree(drSqr_);
	for (int d = 0; d < QmcParams::dim; d++) {
		//ippsFree(resFProposed_[d]);
		//ippsFree(resFOriginal_[d]);
	}
#endif
}
#ifdef POT_LOG
void Interpolate::calcVF() {
	DoublePointerVec resF(QmcParams::dim);
	DoublePointerVec pos(QmcParams::dim);
	dr_.getDimP(pos);
	if (OP_ == ORIGINAL)
	resF = resFOriginal_;
	else
	resF = resFProposed_;
	sizeDimEven_ = sizeDim_ - sizeDimOdd_;
	calcGreen < 3
	> (pos[0] + sizeDimOdd_, pos[1] + sizeDimOdd_, resV_ + sizeDimOdd_, LogArg_
			+ sizeDimOdd_, sizeDimEven_);
	if (numOfOddSlices_ > 0)
	calcGreenForce < 3
	> (pos[0], pos[1], resV_, LogArg_, resF[0], resF[1], sizeDimOdd_);
	for (int m = 0; m < numOfOddSlices_; m++) {
		{
			int zeroPos = m * QmcParams::NumP + zeros_[m];
			for (int d = 0; d < SIM_DIM; d++)
			resF[d][zeroPos] = 0;
		}
	}
	for (int m = 0; m < numOfSlices_; m++) {
		int zeroPos = m * QmcParams::NumP + zeros_[m];
		resV_[zeroPos] = 0;
		LogArg_[zeroPos] = 1;
	}
}
void Interpolate::calcEnergy() {
	sizeDim_ = QmcParams::NumP * numOfSlices_, sizeDimOdd_ = QmcParams::NumP
	* numOfOddSlices_, sizeAll_ = sizeDim_ * SIM_DIM, sizeDimEven_ =
	sizeDim_ - sizeDimOdd_;
	this->calcVF();
}
void Interpolate::calcSingleGreen(BeadPos const &dr, double &v, ForceData &f) {
	double x = dr.getD(0), y = dr.getD(1), logArg;
	calcGreenForce < 3 > (&x, &y, &v, &logArg, &f[0], &f[1], 1);
	v += -0.5 * log(logArg);
}
void Interpolate::calcSingleGreen(BeadPos const &dr, double &v) {
	double x = dr.getD(0), y = dr.getD(1), logArg;
	calcGreen < 3 > (&x, &y, &v, &logArg, 1);
	v += -0.5 * log(logArg);
}
#endif
#ifdef POT_AZIZ
void Interpolate::calcEnergy() {
	sizeDim_ = QmcParams::NumP * numOfSlices_, sizeDimOdd_ = QmcParams::NumP
			* numOfOddSlices_;
	this->calcDist();
	double rZeroFactorT = QmcParams::rZeroFactor;
	double rZeroScaleT = QmcParams::rZeroScale;
	for (int m = 0; m < numOfSlices_; m++) {
		int drSqrPos = m * QmcParams::NumP + zeros_[m];
		drSqr_[drSqrPos] = 10.;
	}
	int index;
	DoublePointerVec pos;
//	if (OP_ == ORIGINAL)
//		std::transform(resFOriginal_.begin(), resFOriginal_.end(), resF.begin(),
//				[](vector<double> &fVec) {return fVec.begin();});
//	else
//		std::transform(resFProposed_.begin(), resFProposed_.end(), resF.begin(),
//				[](vector<double> &fVec) {return fVec.begin();});
	dr_.getDimP(pos);
#pragma omp parallel for
#pragma ivdep
	for (int i = 0; i < sizeDimOdd_; i++) {

#ifdef FORCE
		calcGreenForce(drSqr_[i], rZeroScaleT, rZeroFactorT, resV_[i],
				resF_[i]);
		//}
		for (int d = 0; d < SIM_DIM; d++)
		//double *resFP = resF[d], *posP = pos[d];
#pragma ivdep
		if (OP_ == ORIGINAL)
		//	for (int i = 0; i < sizeDimOdd_; i++)
		resFOriginal_[d][i] = pos[d][i] * resF_[i];
		else
		resFProposed_[d][i] = pos[d][i] * resF_[i];
	}
#else
		calcGreen(drSqr_[i], rZeroScaleT, rZeroFactorT, resV_[i]);
#endif
#pragma omp parallel for
#pragma ivdep
		for (int i = sizeDimOdd_; i < sizeDim_; i++) {

			calcGreen(drSqr_[i], rZeroScaleT, rZeroFactorT, resV_[i]);

		}
		for (int m = 0; m < numOfSlices_; m++) {
			int zeroPos = m * QmcParams::NumP + zeros_[m];
			int drSqrPos = m * QmcParams::NumP + zeros_[m];
			resV_[zeroPos] = 0.;
		}
#ifdef FORCE
		for (int m = 0; m < numOfOddSlices_; m++) {
			{
				int zeroPos = m * QmcParams::NumP + zeros_[m];
				for (int d = 0; d < QmcParams::dim; d++)
				if (OP_ == ORIGINAL)
				//	for (int i = 0; i < sizeDimOdd_; i++)
				resFOriginal_[d][zeroPos] = 0.;
				else
				resFProposed_[d][zeroPos] = 0.;
			}
		}
#endif
	}
void Interpolate::calcDist() {
	dr_.norm2(drSqr_);
}
void Interpolate::calcSingleGreen(BeadPos const &dr, double &v, ForceData &f) {
	double drsqr = dr.norm2(); // / QmcParams::rZeroFactor;
	double fr;

	calcGreenForce(drsqr, QmcParams::rZeroScale, QmcParams::rZeroFactor, v, fr);
	for (int d = 0; d < QmcParams::dim; d++)
	f[d] = dr.getD(d) * fr;//;fgrid_[index];
	//}
}
void Interpolate::calcSingleGreen(BeadPos const &dr, double &v) {
	double drsqr = dr.norm2(); /// QmcParams::rZeroFactor;
	calcGreen(drsqr, QmcParams::rZeroScale, QmcParams::rZeroFactor, v);
}
__declspec (vector (uniform(drTh,rZero))) void calcGreenForce(double dr, double rZero, double &v, double &f) {
	const double D = 1.241314;
	const double A = 0.5448504E6;
	const double EPS = 10.8;
	const double ALPHA = 13.353384;
	const double C[3] = {1.3732412, 0.4253785, 0.1781};
	const double Cd[3] = {1.3732412 * 6, 0.4253785 * 8, 0.1781 * 10};
	double h, p, dp, dh, drsqr, dr_minus_ten, dr_minus_eleven, exp_arg, drS;
	drS = dr / rZero;
//	if (drS > drTh) {
//		v = 0;
//		f = 0;
//	} else {
		dr_minus_ten = pow(drS, -10);
		dr_minus_eleven = dr_minus_ten / drS;
		drsqr = drS * drS;
		p = dr_minus_ten * (C[2] + drsqr * (C[1] + C[0] * drsqr));
		dp = -dr_minus_eleven * (Cd[2] + drsqr * (Cd[1] + Cd[0] * drsqr));
		v = EPS * (A * exp(-ALPHA * drS));
		f = EPS * (-ALPHA * A * exp(-ALPHA * drS));
		if (drS > D) {
			v -= EPS * p;
			f -= EPS * dp;
		} else {
			exp_arg = (D / drS - 1);
			h = exp(-pow(exp_arg, 2));
			dh = 2 * D * h * exp_arg / (drsqr);
			v -= h * EPS * p;
			f -= EPS * (dh * p + h * dp);
		}
		f /= (dr * rZero);
	//}
}
__declspec (vector (uniform(drTh,rZero))) void calcGreen(double dr, double drTh,
		double rZero, double &v) {
	const double D = 1.241314;
	const double A = 0.5448504E6;
	const double EPS = 10.8;
	const double ALPHA = 13.353384;
	const double C[3] = {1.3732412, 0.4253785, 0.1781};
	double h, p, dp, dh, drsqr, dr_minus_ten, exp_arg, drS;
	drS = dr / rZero;
//	if (drS > drTh)
//	v = 0;
//	else {
		dr_minus_ten = pow(drS, -10);
		drsqr = drS * drS;
		p = dr_minus_ten * (C[2] + drsqr * (C[1] + C[0] * drsqr));
		v = EPS * (A * exp(-ALPHA * drS));
		if (drS > D) {
			v -= EPS * p;
		} else {
			exp_arg = (D / drS - 1);
			h = exp(-pow(exp_arg, 2));
			v -= h * EPS * p;
		}
//	}
}
#endif
#ifdef EXT_POT
void Interpolate::calcEnergyExtPot(BeadPos Pos, double &v, ForceData &f) {
	double x = Pos.getD(0), y = Pos.getD(1), sqr;
#ifdef LATTICE
	v = -QmcParams::latticeStrength * 0.5
	* (cos(2 * PI * QmcParams::period * x)
			+ cos(2 * PI * QmcParams::period * y));
	f[0] = QmcParams::latticeStrength * PI * QmcParams::period
	* (sin(2 * PI * QmcParams::period * x));
	f[1] = QmcParams::latticeStrength * PI * QmcParams::period
	* (INV_TAN_T * sin(2 * PI * QmcParams::period * x)
			+ INV_SIN_T * sin(2 * PI * QmcParams::period * y));
#endif
#ifdef STEPS
	int index = this->interpIndex(Pos, QmcParams::gridResExtPot);
	v = vgridExtPot_[index];
	f[0] = f[1] = 0;
#endif
#if not (defined(STEPS) or defined(LATTICE))
	int index = this->interpIndex(Pos);
	v = vgridExtPot_[index];
	f[0] = fgridExtPot_[0][index];
	f[1] = fgridExtPot_[1][index];
#endif
}
void Interpolate::calcEnergyExtPot(BeadPos Pos, double &v) {
	double x = Pos.getD(0), y = Pos.getD(1), sqr;
#ifdef LATTICE
	v = -QmcParams::latticeStrength * 0.5
	* (cos(2 * PI * QmcParams::period * x)
			+ cos(2 * PI * QmcParams::period * y));
#endif
#ifdef STEPS
	int index = this->interpIndex(Pos, QmcParams::gridResExtPot);
	v = vgridExtPot_[index];
#endif
#if not (defined(STEPS) or defined(LATTICE))
	int index = this->interpIndex(Pos);
	v = vgridExtPot_[index];
#endif
}
#endif
